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ABSTRACT 

This paper presents an overview of recent applications of Eulerian-Lagrangian computa- 
tional schemes in simulating transonic flutter instabilities, n this approach, the fluid-structure 
system is treated as a single continuum dynamics problem, by switching from an Eulerian to a 
Lagrangian formulation at the fluid-structure boundary. This computational approach effec- 
tively eliminates the phase integration errors associated with previous methods, where the fluid 
and structure are integrated sequentially using different schemes. The formulation is based on 
Hamilton’s Principle in mixed coordinates, and both finite volume and finite element discreti- 
zation schemes are considered. Results from numerical simulations of transonic flutter insta- 
bilities are presented for isolated wings, thin panels, and turbomachinery blades. The results 
suggest that the method is capable of reproducing the energy exchange between the fluid and 
the structure with significantly less error than existing methods. Localized flutter modes and 
panel flutter modes involving traveling waves can also be simulated effectively with no a priori 
knowledge of the type of instability involved. 


NOMENCLATURE 

a = speed of sound; also location of elastic axis 

c = 2b = airfoil chord 

e - internal energy 

£ = total energy 

/ = body force 

h = bending deflection, positive down 

k = (ob/U= reduced frequency 

= typical section bending stiffness 
K a - typical section torsional stiffness 
m - mass per unit span 

M = Mach number 

p = pressure 
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q, = generalized coordinates 

r = position vector 

r a = nondimensional radius of gyration about EA 
s = blade spacing along leading edge locus 
s, = 5 sin0 

s 2 -S COS0 

t = time 

T = stress vector 

T = kinetic energy 

x„ = nondimensional CG-EA offset 

u,v = velocities in x,y directions 
h = material velocity vector 

U = mesh velocity vector 

U = UJb(O a 

U = strain energy 

£/_ = free-stream velocity at upstream infinity 

a = angle of attack; also torsional deflection 
6 = stagger angle; also node rotation 

p = m Inpb 2 = mass ratio 

p = air density 

o = interblade phase angle 

to = circular frequency, rad/s 

to* = uncoupled frequency in bending 
( 0 a = uncoupled frequency in torsion 

Superscripts and Subscripts 
E = elastic 

F = fluid 

<*> = conditions at upstream infinity 

e = conditions at exit 


INTRODUCTION 

The classical analytical methods and computational schemes for solving fluid-structure 
interaction problems were established in a different era, before the dawn of the computer age. 
Early developments in the field of unsteady aerodynamics, especially the exact solution of 
Theodorsen [1] for two-dimensional incompressible flow over an airfoil executing simple har- 
monic motion, undoubtedly influenced the direction of research. Thus, it not surprising that 
modal methods in the frequency domain played such a central role in the theoretical treatments 
of flutter calculations in the classical and modem texts on aeroelasticity [2-5]. 

Recent advances in supercomputers and massively parallel machines have had a significant 
impact on the feasibility of carrying out nonlinear transonic flutter calculations in the time- 
domain. It was pointed out by the author in Ref. 6 that this increased use of computers to simu- 
late the behavior of physical systems has focused attention on the need to reexamine the exist- 
ing classical approaches to certain problems. The frequency domain is a natural setting for 
problems amenable to the Fourier-type solution methods characteristic of classical unsteady 
aerodynamics, but not for nonlinear aeroelastic problems where the aerodynamic forces are 
computed in the time domain. Although time-domain aeroelastic calculations have by now 
become commonplace, the computational approach followed in most of these studies are along 
classical lines. That is, the fluid and the structure are modeled separately, and then coupled by 
specifying the kinematic boundary conditions at the fluid-structure boundary. The kinetic or 
natural boundary conditions are not treated as boundary conditions in the ordinary sense; 
instead, they provide forcing terms in the governing equations of motion for the structure. 
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This classical approach appears, at first glance, as a natural extension of the analytical 
frequency-domain methods into the computational time-domain. But on closer examination, a 
most important and fundamental difference emerges. Because the boundary between the fluid 
and the structure is a Lagrangian surface, whose state (displacements and velocity) must be 
known in order to impose the kinematic boundary conditions at this surface, the classical time 
marching method cannot be carried out without introducing approximations. The reason why 
approximations are necessary follows from the observation that in order to determine the exact 
state of the boundary at time t, one must first solve the entire system of equations for the struc- 
ture. But this cannot be done until the surface pressure is known, which depends on the solu- 
tion for the fluid domain, which in turn depends on the (yet unknown) boundary conditions 
during the current time step At. Thus, there can never be exact simultaneity between the 
unsteady events calculated in the structural and fluid domains. Since the fluid and the structure 
are in effect integrated sequentially during each time step, often using different numerical 
schemes, undesirable phase and integration errors are introduced. 

There is thus a fundamental difference between aeroelastic simulations, where the structure 
is free to move in accordance with Newton’s Laws, and steady or unsteady flow calculations, 
where the structure is fixed or forced to move in a kinematically constrained manner. This fact 
is seldom mentioned in papers dealing with computational fluid dynamics (CFD), or with the 
computational aspects of aeroelasticity. In unsteady flow calculations where the motion of the 
airfoil is prescribed, the future state of the fluid-structure boundary is known and numerical 
errors (whatever their source) do not affect the future position and motion of the boundary. 
Because the motion of the structure is unaffected by conformity or consistency errors arising at 
the fluid-structure boundary, such errors tend to "average out" over a cycle of oscillation. This 
is not the case in aeroelastic calculations, however. Not only is the future state of the fluid- 
structure boundary unknown, causing the aforementioned difficulty in implementing the boun- 
dary conditions, but errors from all sources affect the future position and velocity of the fluid- 
structure boundary. The computational model of the aeroelastic system behaves as a feedback 
control system, where the state of the structural system is not only affected by the "true" aero- 
dynamic forces, but also by fictitious "control forces" arising from the continued feedback of 
modeling and numerical errors. Here, errors do not tend to "average out" over a cycle and 
remain bounded in the long run, because of the highly nonlinear amplification possible through 
this feedback mechanism. 

Results from a large number of numerical simulations of transonic flutter instabilities indi- 
cate that two distinctly different manifestations of these errors can occur 

1. A long-term cumulative error in the calculated motion of the aeroelastic system, similar to 
the temporal nonuniformities or far-field nonuniformities that occur in perturbation prob- 
lems that are singular or involve multiple time scales. 

2. A significant qualitative and quantitative error in the stability behavior of the system, includ- 
ing a large error in the calculated flutter boundary and possibly also a qualitatively incorrect 
assessment of the type of instability behavior (flutter vs. divergence). 

It should be clear from these observations that one must distinguish between the inherent 
accuracy of the individual CFD and FEM computational schemes used for the fluid and solid 
domains, and the accuracy of the overall aeroelastic scheme. It is meaningless to insist on 
higher-order CFD/FEM schemes in flutter calculations if the boundary conditions are imple- 
mented such that only first-order (or less!) accuracy is obtained in the aeroelastic simulations. 
The implications for implicit schemes, where the allowable time step At is set by accuracy 
rather than numerical stability considerations, are obvious. 

In Ref. 6, we have proposed a new approach to the problem, which has been further 
developed and explored in a series of papers [7-11]. By formulating the governing equations 
for both the fluid and the structure in integral conservation-law form based on the same mixed 
Eulerian-Lagrangian description, the entire fluid-structure system can be treated as one prob- 
lem, while still allowing for different spatial discretization schemes in the two domains, if so 
desired. In this approach, the boundary conditions become transparent to the time-marching 
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algorithm, resulting in highly vectorizable codes. The method also provides a systematic and 
efficient procedure for coupling finite element structural codes to finite volume or finite ele- 
ment fluid codes, allowing the direct use of highly accurate finite element structural models 
rather than a truncated set of normal modes. 

The objective of the present paper is to present a brief overview of the computational 
method, and to illustrate its application in distinctly different areas of classical and computa- 
tional aeroelasticity. The transonic flow regime is intentionally emphasized, because of the 
challenges it represents both from a computational fluid dynamics standpoint, and also because 
of its rich source of nonlinear and surprising aeroelastic phenomena. These cases include 
illustrations of aeroelastic instabilities where several degrees of freedom participate to generate 
localized response or traveling waves, which are difficult to treat accurately using modal 
methods. 


VARIATIONAL FORMULATION OF AEROELASTIC CONSERVATION LAWS 
Hamilton's Principle 

Coupled fluid-structure interaction problems involve, by definition, interactions between two 
distinct material domains (fluid and solid). In modeling the structural problem, finite element 
methods based on a Lagrangian (material) description are typically used, and Hamilton’s Prin- 
ciple and variational methods play an important role in formulating the equations of motion. 
For the fluid problem, finite difference or finite volume methods based on an Eulerian (spatial) 
description are more popular, and variational methods are conspicuous in their absence. 

If Eulerian coordinates are used in connection with the classical Lagrangian density, 
Hamilton’s Principle gives trivial results. Various attempts have been made to augment the 
Lagrangian by adding terms that, when substituted into Hamilton’s Principle will yield the 
sought field equations; see, for example, the discussion in Ref. 12. Although many of these 
investigations take Hamilton’s Principle as a starting point, they usually have little if anything 
to do with the physical contents of the principle. 

Hamilton’s Principle represents an alternate formulation of Newtonian mechanics, and is 
therefore capable of producing the correct equations of motion in an arbitrary coordinate sys- 
tem. There is no need to introduce ad hoc terms in the Lagrangian, a procedure that this author 
considers in violation of the spirit of Hamilton’s Principle. One must recognize, however, that 
Hamilton’s Principle is based on material particles and variations of the Newtonian paths of 
such particles, and that the mass of a particle is not to be varied in the process of varying the 
path. 

To this end, we introduced in Ref. 7 the concept of a "material variation”, constructed such 
that the variational operator 6 in Hamilton’s Principle remains a Newtonian path variation 
under coordinate transformations. Consider a finite region of a continuum, represented by a 
materia] cell of volume Q, as illustrated in Fig. 1. Hamilton’s Principle in its extended form, 


*2 

J [ST-W +&W E ]dt=0 


( 1 ) 


is applicable without additional assumptions regarding the nature of the forces present. Here, 
5W e is the virtual work of all forces acting on the material particles in the volume, and may 
include forces for which no scalar potential can be found such that SW E = -5P. In this very 
general form, the principle is an integral statement of D’Alembert’s Principle rather than a 
"true" variational principle. 

If all forces are derivable from scalar potentials, Eq. (1) becomes a variational principle: 
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( 2 ) 


ti n 

b]Ldt = bjj L 0 dQdt = 0 

<i t\Q 

This form is well-established for conservative dynamical systems; however, it is also valid for 
nonconservative systems where the Lagrangian is time-dependent, as long as the previous 
assumptions hold. 

The Lagrangian density L 0 is defined as 

l 0 = t 0 -u 0 -v 0e = t 0 -Vq ( 3 ) 

where 

T 0 = J pTiiiii + (4a) 

U 0 =jo ij de ij (4b) 

0 

Here, T 0 and U 0 are the kinetic and internal energy densities, respectively. The terms 7 q and 
To are only present in noninertial reference frames and give rise to centripetal and Coriolis 
forces. The coordinates t|, are arbitrary generalized displacement coordinates, in the distri- 
buted sense. Finally, V 0e is the density associated with the potential of the external forces, 

VVf = J VV 0e da (5a) 

a 

p/=-W 0£ (5b) 

For a solid element, U o is the familiar strain energy density. For a fluid element, U o is equal to 
the (intrinsic) internal energy per unit volume, which for a perfect gas reduces to 

(6a) 
(6b) 

where e is the total energy per unit mass. In either case, U o represents the internal potential 
energy, and Eq. (6a) is obtained from Eq. (4b) by assuming inviscid flow, = -p5y, and 
integrating along a reversible path, using the perfect gas law. 

Euler-Lagrange Equations of Motion 

For inviscid flows, it suffices to consider variations of "action" integrals of the form 


U „ = 


Y-l 


pe = t/o + Tq 


<2 


/= J j L 0 (x i ,t,i\ h r\ iJ ,i\i)dadt 

• l 0(0 


(9) 


where the ti.’s represent arbitrary generalized displacement coordinates in the distributed 
sense. In Lagrangian coordinates, applying Hamilton’s Principle by setting 8/ = 0 ( £i fixed) 
leads to the classical Euler-Lagrange equations of the calculus of variations: 
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( 10 ) 


3L 0 3 3Lq ^ ^ dLp ^ _ q 

3th 3/ dt|. 3 j tj dr\ij 


In a mixed Eulerian-Lagrangian description, the domain Q(f) constitutes an intermediate refer- 
ence frame, moving with an arbitrary but kinematically specified velocity U . Consider the 
material integral 


(HD 


(ID 


The time rate of change of <I> in the jc, reference frame is 

d<t> 


p _ (• 

dt (HD 


& + V-[(« - 1/>I>] 
ot 


(12) 


and may be considered a total or convective rate. It is then straight-forward to show that, in 
order to preserve the material nature of the variational operator ”8" in Hamilton’s Principle, the 
time derivative in the classical Euler-Lagrange equations must be changed as follows: 


^(•)-»|’0 + v [( «-i / X ) ] 


(13) 


The corresponding (modified) Euler-Lagrange equations, valid in arbitrary mixed cartesian 
coordinates, become 


It - J- ^ - T- |( “;- y ^ + I t 1 " 0 

dn. dl dni Bx j ^ ®nij 


(14) 


In an inertial coordinate system x, in the absence of body forces, dL 0 ldr\i = 0 (translational 
invariance). The following conservation-law form of Eq. (14) is then obtained: 


3 df'O 

* dri7 + 





dLp 

dn.j 


] = o 


Recognizing that the velocity field is given by «, = 


dL o 


= P/i 



and because of the symmetry of the strain tensor £,y, it follows that 


dL 0 
dn ij 


dUp Wp 


(15) 


(16) 


(17) 


The governing equations (14) can thus be written as 
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(18) 


(puj) + Vj)pi*i ~o,j)] = Pfi 

valid in an inertial cartesian coordinate system x ,. For an inviscid fluid, the stress tensor is 
diagonal 

<*ij = ~P&ij < 19 > 

and we obtain the Euler equations of fluid mechanics by setting Uj = 0. If one sets u ; = Uj, the 
Lagrangian form of the equations of motion are obtained (the divergence is readily expressed 
in terms of Lagrangian coordinates; see, for example. Lamb [13].) 

We note in passing that the corresponding Navier-Stokes equations can be obtained by sub- 
stituting into Eq. (18) the constitutive law 


o,j = -p% + 2p( e.j - y A8y) 

(20) 

where p here is the viscosity and 


&ij — ~2 ^ u U + ^ 

(21) 


are the strain rates and the volumetric expansion rate, respectively. 

The conservation of mass equation is a kinematic constraint, and is readily incorporated into 
the variational formalism using a Lagrange multiplier. In the present paper, the integral form 
is used in the actual discretization procedure, 

4- f pd(l + J p («; - Ui) n, dS = 0 (22) 

™ h an 

In most cases in aeroelasticity, the energy equation can be handled in a similar manner, since 
in cases where only mechanical energies are involved, it is simply a first integral of the equa- 
tions of motion. In integral form, it reads 

4-\pedCl+ l pe (u, - U t ) n t dS = \ puj, dQ. + f u,77 dS (23) 

ot ] Q & a a o 

where / and 77 are the cartesian components of the body forces and surface tractions, respec- 
tively. The components 77 are related to the stress tensor through the relation 


77 — Ojy flj 

It is convenient and useful to integrate the energy equation and use the relation 


(24) 


p = (Y-l)[pe-T 0 ] (25) 

to eliminate the pressure p. For a solid cell, such an elimination is not generally possible, 
because the internal energy depends on all six components of the stress tensor. 
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Boundary Conditions 

In the c lassi cal method of numerically solving fluid-structure interaction problems, only the 
kinematic boundary condition of tangent flow is imposed: 


+ u-VB = 0 (26) 

dt 


where B(x,y,z,t) = 0 defines the instantaneous locus of the fluid-structure boundary. The 
kinematic boundary conditions are not enforced explicitly; that is, one does not enforce local 
force equilibrium between a fluid and a solid element at the boundary B, nor does one equate 
the energy flow or power at the common boundary . 

In the present study, both kinematic and kinetic boundary conditions are satisfied locally at 
the fluid-structure boundary. By switching from an Eulerian to a Lagrangian formulation at 
the boundary, the boundary conditions become transparent and they can be satisfied automati- 
cally as the equations of motion are integrated. The kinematic boundary condition can be 
imposed by zeroing out the convective fluxes at the fluid-structure boundary, leaving only the 
pressure terms in the boundary integrals. These pressure terms contribute to the kinetic boun- 
dary conditions, but need not be considered explicitly because the kinetic boundary conditions 
between cells are taken care of by the conservation laws being integrated. 

At the upstream and downstream far-field boundaries we use nonreflecting boundary condi- 
tions of the type formulated by Hedstrom [14] and generalized to two space dimensions in Ref. 
15. In the cascade calculations, phase-lagged periodic boundary conditions are applied at the 
boundaries of the reference channel. Thus, if % denotes the spatial coordinate along the 
periodic (upper and lower) boundaries of the reference channel, we enforce 


W l a,t)=W u a,t-ola) (27a) 

(27b) 

where a and to are the interblade phase angle and the circular frequency, and the subscripts u 
and / indicate upper and lower boundaries, respectively. In order to prevent numerical instabil- 
ities, it is necessary to limit the magnitude of the change in the periodic boundary conditions 
during each time step. 


SPATIAL DISCRETIZATION SCHEMES 
Finite Volume Method 

In the two-dimensional cases considered in this paper, we take the unknowns as 



P 

P“ . 

pv 

pe 


(28) 


Applying the divergence theorem to the Euler-Lagrange equations, with the mass and energy 
equations appended, the system is transformed into the following set of integral conservation 
laws: 
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(29) 


d_ 

dt 



n x + Fytly )ds = 0 


where ft is an element area with (moving) boundaries dft and 



P(u-U) 


p(v-VO 


pu(u-U)-a„ 

L r __ „ 

p«(v-V) + a„ 

F x =‘ 

pv(u-U)-a„ 


pv(v-V)-a„ 


pe(u-U)-a„u-a„v 


pe(y-V) + a„u-a„v 


Here, u,v and U.V are the cartesian components of u and U, respectively, and a„ and o„ are 
the normal and shear stresses at the surface 3ft. For a cell occupied by an inviscid fluid, we set 
= ~P> = 0 and use Eq. (25) to eliminate the pressure. 

By applying the integral form of the conservation laws to each fluid cell, we obtain a system 
of coupled ordinary differential equations of the form 

+ Qij - D 0 = 0 (31) 


where fty is the cell area, W tJ is the vector of unknowns for cell (ij), and Qy is the net flux 
out of the cell (i,/), contributed by the integral over 3ft in Eq. (29). The flow variables are 
assumed constant over each cell and are averaged at cell edges, resulting in second-order accu- 
racy on a smooth mesh. D tJ is a dissipative operator added to damp out numerical oscillations 
and to prevent decoupling of even and odd cells. The dissipative fluxes are constructed accord- 
ing to the idea of "adaptive dissipation" developed by Jameson and Baker [16], using a blend 
of second and fourth differences in the flow variables. 

Galerkin Finite Element Method 

Two different finite element schemes are considered in this paper. The first scheme is closest 
to the finite volume procedure discussed in the previous section, wherein the divergence 
theorem is used (in reverse) to obtain the integral conservation-law form of the equations. The 
Galerkin finite element procedure is then applied to obtain space-discretized equations in the 
usual manner, by approximating the vector W of unknowns as 


W(x,yj) = ^Wj(t)^(xj) (32) 

j 

where Wj are nodal values of W and <[>, are shape (interpolation) functions. In the present 
paper, linear interpolation functions are used and the solution is implemented using triangular 
elements. 

The solution is obtained by discretizing the weak form of the equations, 

-|-JW'(Mft + !<t> i V-Fdft = 0 (33) 

where 

F = F x i+F y j (34) 
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and F x and F y are given by Eq. (30). Substituting Eq. (32) into the first term in Eq. (33) and 
interchanging the order of summation and integration, one obtains 


where 

m% = l$Ajd £2 


(35) 

(36) 


are the elements of the consistent mass matrix, for element £2. In the unsteady (deforming 
mesh) case , the mass matrix is time-varying and depends on the local Jacobian determinant 
However, if linear shape functions are used, the generalized mass elements can be evaluated 
exactly. 

Using the divergence theorem, the second term in Eq. (33) can be written as 

kV F<JQ= f hF ndS-jF Vhdn (37) 

a da a 

where n is the unit outward normal to the boundary 3£2 of £2. Because F is a nonlinear func- 
tion of the flow variables, the integrals in Eq. (37) must be evaluated numerically. Denoting 
the contribution by the flux integrals, Eq. (37), associated with node i of element £1 by CP » the 
space-discretized FE equations of motion for a single element become 

j-(? t m l ijWj) + Q? = 0 (38) 

at j 

where the range of the i and j indices are the nodes for element £2. 

Equations (38) are formally identical to the corresponding FE equations for a structural ele- 
ment, if the Wys and the QP’s are interpreted as the element generalized coordinates and the 
corresponding (nodal) generalized forces, respectively. But in the present formulation, we do 
not assemble the element matrices into global system matrices, as is customary in structural 
dynamics. Instead, we perform a local assembly only, by algebraically summing the contribu- 
tions from all elements £2* that meet at node i. The union of these elements forms a control 
volume V, surrounding node i; hence this local assembly can also be performed by applying 
the above finite element discretization equations directly to the "superelements" 

V,= I£2i ; i-1, 2, • • ■ ,N (39) 

*=i 

where Ni is the number of elements meeting at node i and N is the total number of nodes. 
Adding dissipative terms, the space-discretized Galerkin FE equations can be written as 

j-('Lm ij Wj) + Q l -D i = 0 (40) 

at j 

where i = 1,2, • • • N (all nodes), and the summation on j extends over the nodes in the 
superelement V i associated with node i. The dissipative fluxes D, are of the Jameson type, and 
are slight modifications of the terms used by Mavriplis [17] in his nodal finite volume code. In 
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the FE calculations presented in this paper, unstructures triangular grids of the type illustrated 
in Fig. 2 were used. Slightly different implementations were used in the external (panel) and 
the internal (cascade) calculations; see Refs. 18-19 for further details. 

Lagrangian Finite Element Method 

The second scheme proceeds directly from the variational principle (Hamilton’s) to the 
discretized Lagrange equations in a suitable set of generalized coordinates, (<?)y. This is the 
procedure used to obtain the space-discretized equations for the solid domain (wing, blades, 
panel). At the fluid-structure boundary, we switch from a mixed method to a Lagrangian 
method by setting C/* = u<. 

Modified Typical Section. A typical section isolated wing model has been implemented as 
illustrated in Fig. 3. A section of unit width in the spanwise direction is considered, in the 
spirit of the original ideas of Theodorsen and Garrick, but with two important differences. 
First, the section is allowed to have camber bending, and this "chordwise" flexibility is 
modeled using plate-type elements of unit width; see Fig. 3b. Each element is allowed a hol- 
low core, by specifying the effective structural skin thickness t, for each element. Second, 
since the method of calculation is at the element level, an attempt is made to model the distri- 
buted restraining stiffness from the remainder of the wing on each element, by introducing 
bending and torsion springs at the element nodes. By careful tailoring of the parameters as dis- 
cussed in Ref. 6, the strip model can be matched to any regular typical section model in the 
limit as the chordwise stiffness approaches infinity. 

On a C-mesh of quadrilateral elements, one obtains a set of discretized equations very simi- 
lar to Eqs. (31): 


^((m]y{9) iy )+ {Q E hj-{Q F hj = 0 (41) 

Here, the ij subscripts refer to the (ij) node, whereas elements are labeled with superscripts; 
i.e., [m]** is the element mass matrix. For the blade structure in Fig. 3, the structural elements 
are at ;=1 , and we can drop the j - subscript and write 

Wlf={". 9. ) (42) 


where w>, and 0, are the transverse and rotational displacements at node (i, 1), Fig. 3b. Plate 
element of unit width are used, with the transverse displacement w(£,/) inside a typical ele- 
ment approximated in terms of the nodal coordinates q k as follows: 

w&O-IAWBfcO 0 ( 43 ) 

*=i 


where the N k ’s are the standard cubic shape functions. 

The 2x2 "nodal" mass matrix in Eq. (41) is written in terms of the lumped element mass 
matrices for elements i and i+1 as 


[m], = [mYn + [mJn 1 


(44) 


where the element matrices have been partitioned into 2x2 submatrices in the usual way. 
Similarly, the nodal elastic forces (G £ )i>. and die generalized nodal fluid pressure forces 
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{Q F )ij can be expressed in terms of the corresponding element forces, as follows: 


where 


Gw=Gi <i+1) + 6? 

Qi=Q F 2 iM) +Q F A i 




(45) 

(46a) 

(46b) 

(47) 


and p) and pi are the fluid pressures on the lower and upper surfaces of the ith element, respec- 
tively. It is important to note that we do not need to assemble the structure and form global 
stiffness (or mass) matrices; only local assembly at the element (node) level is required, pre- 
cisely as in the fluid domain. Although we have assumed a linearly elastic structure in Eq. 
(45), nonlinear material behavior, both elastic and inelastic, are relatively easy to incorporate 
at the local element level. 

To complete the solution procedure, we need to perform another integration to get from {q } 
to [q ) . This is necessary since we need to know the actual nodal displacement vectors in order 
to move the mesh coordinates in a Lagrangian fashion. It is done simultaneously in the 
multi-stage Runge-Kutta scheme, by writing the equations in state-variable form by appending 
at each node the matrix identity 

«*> 

The energy equation for the structure is discretized and integrated as follows 

(49) 

dt 


where 

+ (50) 
i 

For a linearly elastic structure, the kinetic and strain energies can be expressed in terms of the 
present (local) element mass and stiffness matrices as 

Ti + Ui = j[q )f[m]i tfh + H *],{?), (51) 

From a computational standpoint, it is both natural and advantageous to use finite element 
throughout both the structural and the fluid domains. By imposing "conformity" or "compati- 
bility" conditions on the shape functions, we guarantee that the basic conservation laws for 
mass, momentum, and energy are satisfied at the fluid-structure boundary, to within the discret- 
ization error of the scheme. For maximum modeling flexibility and computational efficiency, 
both the fluid and the solid element should be allowed interior nodes with respect to the other. 
In modeling thin compressor blades with plate elements, as in the present scheme, we use a 
much coarser mesh for the structure, especially in the leading and trailing edge regions where 
flow gradients are high. By matching the mesh densities appropriately, a simultaneous optimi- 
zation of the computational efficiency and the discretization/modeling errors can be achieved. 
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Rotating Turbomachinery Blade. The blade is modeled as a thin plate with varying thick- 
ness, using the von Karin an theory to incorporate the effect of the in-plane centrifugal forces. 
The kinetic energy of the rotating blade. Fig. 4, can be expressed in the form 

r=yJp/» {w 2 + ft 2 [(r 0 +y) 2 + (wcosS + xsinO) 2 ]} ditty (52) 

where term associated with Coriolis forces have been neglected. The total strain energy can be 
expressed as 



U = J U o dxdydz = U b + U a 

(53) 


fM\hr + 2vf? 

[ax 2 J [dy 2 J dx 2 dy 2 L^J . 

(54a) 


U a = ± J [N x t zo +N y ty 0 + N^Jdxdy 

(54b) 

where D = £/i 3 /12(1 - v 2 ) is the flexural rigidity of the plate. U b is the bending strain energy, 
and U a is the strain energy due to the in-plane centrifugal forces. 

The stress resultants N x ,N y , and N v , which arise from centrifugal forces and give rise to the 
strains e, 0 , ty 0 and y^ 0 in the middle surface of the plate, are approximated as follows: 


b 

N x = J p/iQ 2 £ sin 2 0dlj 

X 

(55a) 


L 

N y = j p/i£2 2 (r 0 + r\)di\ 

y 

(55b) 


N xy = 0 

(55c) 

The virtual work done by the aerodynamic forces on the rotating blade is given by 



bW E = j(pi - p u )hwdxdy = 

(56) 


where p, and p u are the fluid pressures on the lower surface and the upper surface of the blade, 
respectively, and Qf are the generalized fluid forces. 

In the present study, the blade is discretized into nxm rectangular plate elements. For each 
element, a bi-cubic polynomial is used to describe the deflection field. The transverse dis- 
placement inside element ij can be expressed in local element coordinates (£,t|) as: 

w*(tV)=XWU)?*(0 < 57 > 

*=i 
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where q x -q^ q*-qn and <7 13 -fl 16 arc the nodal value of w , dwldx , 3 2 w/9x3> 

associated with nodes 1, 2, 3 and 4. The shape functions Af*(£,T|) are products of the first-order 
one-dimensional Hermitian interpolation functions. 

The space-discretized equations of motion for a given element or the entire blade can be 
obtained using Lagrange* s equations, after substituting Eq. (57) into the appropriate expres- 
sions for T and U: 


^-([M 2 }[q)) + (XK}-[M 0 ]){q) - [Q c )~ [Q F ) = Q (58) 

at 

where (?) is the generalized coordinate vector of the discretized system, [M 2 ] is the mass 
matrix associated with the part of the kinetic energy involving the quadratic terms of general- 
ized speeds, while [M 0 ] and {Q c ) are related to the kinetic energy due to the centrifugal 
force. [If ] is the stiffness matrix including both the linear stiffness and nonlinear stiffness due 
to the in-plane stretching. For further details, see Ref. 20. 

The complete aeroelastic model is shown schematically in Fig. 4b. The cantilever blade lies 
in x-y plane, where x indicates the chordwise direction and y indicates the s panwise direction; 
see also Fig. 4a. The shaded strips are the control channels, where we assume that the motion 
of the fluid is governed by the two-dimensional Euler equations in the x-z plane. The aero- 
dynamic loads are computed inside each channel and interpolated linearly between adjacent 
channels. This channel theory can be interpreted as an extented (numerical) version of classi- 
cal strip theory. 

Nonlinear Panel. The structural model used in the two-dimensional panel flutter simulations 
is based on the nonlinear von Karman plate equations, which can be written as 


d ^- <w » +w -)0 + p*I 


(59) 


Here, N x0 is the initial in-plane loading due to external forces, and N x is the additional in-plane 
force induced by the transverse deflections w: 




2c rax 


(60) 


where c is the plate chord length. 

These equations can be discretized using standard finite element techniques, resulting in a 
set of governing equations of the same basic form as Eq. (41). In the present case, the elastic 
forces in the ith element can be expressed as a sum of two terms 


{Q £ M[*] i +[*,] , ){<7}' 


(61) 


where [A] is the linear part of the stiffness matrix and [k,] is the so-called geometric stiffness 
matrix arising from the nonlinear stiffness caused by the in-plane stretching due to transverse 
deflections of the plate. It should be noted that, because of the fact that the geometric stiffness 
matrix is proportional to (N x0 + N x ), this matrix depends on the deformations of all elements in 
the panel, as is obvious from Eq. (60). Expressions for these matrices can be found in Ref. 21. 
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NUMERICAL EXAMPLES „ . ^ . 

Sub-transonic flows (Merit < < 1) are characterized by local regions of embedded super- 

sonic flow on the upper and/or lower surfaces of the airfoil. At the aft end of these regions, the 
flow is decelerated to subsonic speeds through nearly normal shock waves. During vibration 
and flutter, the shocks move along the airfoil surfaces, changing in strength and possibly van- 
ishing over part of the cycle. 

Because of the important role that shock motion plays in the flutter problem, a realistic 
simulation of transonic flutter instabilities cannot be achieved unless the moving shocks are 
treated correctly and in a time-accurate manner. Furthermore, in order to model a mixed 
subsonic-supersonic flow-field with finite-amplitude shock motion, nonlinear field equations 
are required. If the flow is not of the mixed type, i.e., is either entirely subsonic or entirely 
supersonic, the equations can usually be linearized for sufficiently small amplitudes of motion, 
or for sufficiently high reduced frequencies. In the mixed subsonic-supersonic (transonic) 
case, however, the linearization is subject to severe restrictions on blade amplitudes and may 
in some cases break down. 

Transonic Flutter of Airfoils 

Except as noted, we have used 5 finite elements in the chordwise direction; 3 elastic ele- 
ments and two small rigid-body elements at the leading and trailing edges. An all aluminum 
structure was assumed, with t,l t varied as noted on the figures. 

Figure 5 compares the calculated bifurcation diagrams for the NACA 64A006 wing model 
(p = 10), as calc ulate d using the present method and compared with results taken from Fig. 3 
of Ref. 22, which were calculated using the classical approach, as implemented in Refs. 22-23. 
This model has also been studied extensively by Ashley [24]. In the notation used here, it has 
the following nondimensional parameters: 

a = - 0.2 

x a = 0.2 (Xc = °) 

r„ = 0.29 (ph = 0.25) 

oycDa = 0.3434 (lOw/CDe = >/5T) 

The corresponding equivalent parameter values in Ashley’s notation have been indicated in 
parenthesis. 

To provide meaningful comparisons, the fluid domain was modeled and solved in exactly 
the same manner as in Refs. 22-23, using the same mesh and the the same five-stage Runge- 
Kutta integration scheme. However, in the earlier method the structural equations were 
integrated separately, after the fluid equations had been advanced to the next stage within the 
multistage scheme. This necessitated an approximation for the lift and moment for the next 
time step, and a linear extrapolation was used in Refs. 22-23. In the present approach, no such 
approximation is necessary, since the entire domain of the continuum and all its cells are 
advanced simultaneously in time. That is, no distinction is made between solid and fluid cells 
during the actual time integration process. 

From Fig. 5 it is evident that the linear flutter boundary, represented by the bifurcation 
points on the U axis, are in reasonable agreement. However, for this example at least, the clas- 
sical method underpredicts the limit cycle amplitudes for reduced velocities significantly 
beyond the linear flutter boundary. The reason for this is that the method underpredicts the net 
energy flux from the fluid to the structure. Close to the linear flutter points, the results are less 
definite and show a greater sensitivity to the computational mesh. At a Mach number of .87, 
decreasing the chordwise stiffness of the model by decreasing the structural thickness ratio t,H 
has the effect of increasing the energy flow to the structure, resulting in an increase in the 
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flutter limit cycle amplitude. A comparison of the long-time limit cycle amplitudes is shown in 
Fig. 6. This camber bending effect appears negligible for ratios of 0.25 and higher, but 
become significant for t,!t less than 0.1. This is in contrast to the same wing at Af = 0.92, 
where chordwise flexibility initially (close to the flutter boundary) is stabilizing, then becoming 
mildly destabilizing as U is increased beyond 2.3. 

For Mach numbers in the upper transonic range (0.9- 1.0), the shocks are typically found to 
be confined to the trailing edge region of the airfoil. During flutter, the shocks do not move but 
change in strength, imposing a significant localized loading on the trailing edge. If the trailing 
edge is sufficiently flexible, the present model predicts that a local "trailing-edge buzz" mode 
can get excited, as illustrated in Fig. 7, where the trailing-edge region oscillates somewhat like 
a trailing-edge flap. Here, the Mach number is 0.92, U = 3.0, and t,U = 0.10. 

The term "buzz" is used here to suggest a high-frequency modulation of the flutter mode, 
and not necessarily to suggest any connection to the well-known phenomenon of "aileron 
buzz", where the viscosity is generally believed to play an important role, through a shock- 
induced separation of the boundary layer. However, it seems also plausible that the same 
shock-dominated mechanism, albeit inviscid, that here is observed to cause trailing-edge buzz, 
could also produce an aileron buzz phenomenon. This question was addressed in Ref. 11, 
where it was found that such nonclassical (inviscid) buzz phenomena are indeed predicted by 
the present Euler-based calculations, at Mach numbers very close to those where buzz has 
been observed in wind tunnel and flight tests. Figure 8 illustrates such a nonclassical aileron 
buzz instability, as simulated by the present code. 

Figure 9 shows the results of flutter simulations for an NACA 0012 model, with typical sec- 
tion parameters similar to those of the NACA 0012 Benchmark Model tested at NASA Lang- 
ley [25]: 


a = 0.0 


x a = 0.0 


= 0.25 

(0/,/cQq = 0.6564 

Limit cycle flutter is predicted, in a predominantly bending mode, which is in good agreement 
with the observed flutter at this flutter point. Note the very high mass ratio (4284). 

At more realistic mass ratios, our simulations predict a nonlinear interaction between flutter 
and a weak divergence instability at this same Mach number (0.80), as discussed in Ref. 26. 
The basic mechanism behind the emergence of this weak divergence, which is quenched by 
the flutter instability (see Fig. 4 of Ref. 26), is the strong sensitivity of the shocks and the 
supersonic pockets to relatively small airfoil displacements, as illustrated in Fig. 10. This leads 
to a symmetry-breaking bifurcation, where the op=0 equilibrium position becomes unstable and 
bifurcates into two stable nonzero (but small a) equilibrium points. 

An example of the Type 2 error discussed in the Introduction is presented in Figs. 1 1 and 12. 
At lower transonic Mach numbers, reasonable agreement is observed between the calculated 
flutter behavior of the NACA 0012 model, using the classical vs. the present calculation 
approach. In the upper transonic range ( M > 0.9), however, significant discrepancies were 
observed. At Mach 0.94, for example, the present simulations indicate that the linear flutter 
boundary is at around U = 6.1, with strong flutter already evident at U = 6.5, as shown in Fig. 
11. The classical calculation method predicts no flutter whatsoever, for all reduced velocities 
up through 8.0 (see Fig. 12), and predicts divergence around U = 8.2-8.5. A doubling of the 
mesh density and a reduction in the time step used had essentially no effect on the predicted 
behaviors. The precise reason for this significant difference in predicted aeroelastic behaviors 
is presently unknown 
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Transonic Flutter In Cascades and Rotors 

Fan and compressor blades are thin, flexible structures, usually of solid titanium or steel, 
where camber bending is of considerable importance because of the presence of plate-type 
modes. In addition to being susceptible to a number of different flutter instabilities, all tur- 
bomachinery blades are subject to fatigue. It is therefore highly desirable to have flutter and 
aeroelastic response computational methods based on realistic finite element models for the 
blades, so that dynamic stresses can be calculated. 

In the present flutter simulations, two different structural models are used. In the 2D cas- 
cade calc ulatio ns, we make use of the modified typical section model discussed previously. In 
this model, the blades are treated as thin isotropic plates of unit width in the spanwise direc- 
tion, with the typical section stiffness represented by attaching blade bending and torsion 
springs K k and K a at the elastic axis of the blade. The chordwise flexibility (camber bending) 
of the blade is modeled using several plate elements in the chordwise direction, and the fluid 
domain is discretized using the finite volume procedure. In the quasi-3D calculations, we 
make use of the rotating blade model discussed above, and the fluid domain is also solved 
using finite element discretizations. 

Because the blade sections are flexible, the typical section stiffness parameters have been 
nondimensionalized as follows: 


K k = 


12 A 

Et 3 


12 cK a 



(62) 


where t is the maximum thickness of the blade section. Standard sea level values for the inlet 
static pressure and density are used, and the blade material is titanium. Using Eq. (62), the 
bending-torsion frequency ratio can be expressed as 



where r a is the radius of gyration of the blade section about the elastic axis, nondimensional- 
ized by the semichord b. Furthermore, 



3 jcpy 



(64) 


where is the free-stream speed of sound, y is the ratio of specific heats, E is Young’s 
modulus, and p = m /(npb 2 ) is the mass ratio of the blade. 

Unstaggered cascades represent perhaps the best test cases for code checkout and validation 
purposes, because of their simple geometry and boundary conditions. Linear theory predicts 
that the critical interphase angle for flutter of an unstaggered cascade is a = 180 deg; that is, 
adjacent blades are in anti-phase motion with respect to each other. The boundary conditions 
at the channel side boundaries thus become simply the condition for tangent flow, i.e., v = 0. 
This boundary condition remains exact in the nonlinear case, under the assumption that the 
180-deg mode is realizable and corresponds to the least stable interblade phase angle. 

An illustration of limit cycle flutter in cascades is presented in Fig. 13. After a slight 
overshoot, the blade settles into a stable limit cycle. Part a) of this figure shows the transverse 
bending displacement w = -hlc (nondimensionalized by the chord c = 2b) at the elastic axis, 
and the local rotational displacement 0 about this axis. Note that 0 = -a if camber bending is 
neglected. Part b) shows the corresponding nondimensional total energy, E u , = U + T, of the 
blade section vs. time, and the work W A done by the fluid pressure on all the finite elements 
making up the blade section. Also plotted is the difference E ul - IV A , which in the present 
undamped model should equal to a constant, i.e., the initial energy of the blade. By integrating 
the energy equation, we obtain an independent check on the accuracy of our method and the 
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Runge-KuUa integration scheme. Note the excellent agreement between the initial energy and 
the difference E u , - W A , over more than 20 flutter cycles of oscillation, representing roughly 
50,000 time steps. This attests to one of the important advantages of the present method in 
aeroelastic calculations! namely, the fidelity by which the transfer of energy (and therefore also 
momentum) between the fluid and the structure is reproduced. The C p behavior during flutter 
is shown in Fig. 14, which reveals a significant shock motion, of the order of the blade chord. 

Figure 15 illustrates the emergence of a "nonclassical" flutter mode in a staggered cascade 
of NACA 0006 blades, operating at an angle of attack of 3 degrees and at a Mach number of 
0.85. The interblade phase angle is zero. In this example, the bending and torsion modes 
interact nonlinearly to yield a flutter mode that is not a simple, exponentially growing sinusoid. 
Instead, the flutter time history reveals a regular pattern that repeats (with increasing ampli- 
tudes) every 3 cycles. Within the 3-cycle pattern, the bending-torsion amplitude ratio changes 
over a significant range, hlba = 02-0.1. This suggests a nonlinear interaction between the 
bending and torsion modes, wherein energy is continually exchanged between the two modes. 
Although the flutter mode is not a "classical" flutter mode, it is nevertheless a clearly defined 
coherent aeroelastic mode. The blades undergo significant camber bending during flutter. 

For the stagg ered NACA 0006 cascade considered in Fig. 15, an increase in the Mach 
number from 0.85 to 0.90 turned out to be stabilizing, leaving the blade essentially on the 
flutter boundary, as illustrated in Fig. 16. Again, a highly nonlinear flutter mode emerges. 

Finally, an unstaggered cascade of rotating low-aspect-ratio NACA 0003 blades is shown in 
Fig. 17. Such blade sections may be considered representative of the aft stages of compressor 
rotors, or of the thin tip region of advanced fan blades, outboard of the shrouds. Figures 17a) 
and 17b) show the aeroelastic response in bending and torsion, respectively. Here, the angle of 
twist of the blade is defined as 


6 * = ( w te ~ w le) /c 


( 65 ) 


and is equivalent to the conventional torsion theory definition when the blade is chordwise 
rigid. Both figures reveal an increase of the amplitude of motion vs. time, i.e„ flutter, and a 
coherent aeroelastic mode emerges as the blade becomes unstable. The energy plot is illus- 
trated in Fig 17c) and shows that a large amount of energy is extracted by the blade from the 
fluid, indicating explosive flutter. 

Transonic Panel Flutter 

Nonlinear panel flutter in the supersonic range has been studied by a number of researchers, 
using a variety of unsteady aerodynamic theories (for an excellent review of the subject, see 
Dowell [27]). It is well-known that linear aerodynamic theories give reasonable agreement 
with experiments for Mach numbers above about 1.3- 1.4, but give poor agreement in the tran- 
sonic range, 1 < M < 1.3. 

One obvious reason for this discrepancy could be that nonlinear aerodynamic effects are 
important and need be included. Recent simulations using the present computational model 
[10] suggest that the nonlinear aerodynamic characteristics of transonic flows, including mov- 
ing shocks and embedded subsonic regions, may also give rise to a change in the flutter mode 
from a standing to a traveling wave. It was noted by Dowell [27] that in the low supersonic 
(transonic) region, the flutter mode predicted by linearized aerodynamic theories is essentially 
a single-degree-of-freedom mode. However, the present nonlinear calculations using Euler- 
based aerodynamics predict traveling wave flutter modes in the transonic Mach number range. 

Figures 18 and 19 show the panel deflections at 1/4-chord, midchord, and 3/4-chord during 
limit cycle flutter of a typical panel at Mach numbers of 1.0 and 1.2, respectively. The motion 
of the panel at Mach 1 is especially irregular (but periodic) and shows significant phase differ- 
ences between the deflections of the panel at different chord locations, indicating the presence 
of a traveling wave (in the generalized sense used by Dowell). The direction of propagation of 


O. Bendiksen 18 



the wave is in the direction of flow, in agreement with experimental observations. A snapshot 
of the panel deflection and the associated Mach number distribution, Fig. 20, clearly shows the 
presence of a strong shock on the panel. 

When the Mach number was increased to 1.2, the strong shocks disappeared and the limit 
cycle flutter became more "regular”, Fig. 19. However, the presence of higher harmonics in the 
flutter response is still evident, as is the temporal phase shifts between the various points along 
the plate chord. Figure 21 presents a series of snapshot plots of the panel deflection during one 
period of limit cycle oscillation. It is interesting to note that the instantaneous mode shape 
continuously evolves during each cycle of oscillation, in a nonharmonic manner. Again, the 
flutter mode must be classified as a traveling wave, in the generalized sense. 

SUMMARY AND CONCLUSIONS 

An overview of a new method for simulating the aeroelastic response of a diverse class of 
systems has been presented, and the advantages of the new scheme in actual flutter calcula- 
tions have been explored. The method has number of important advantages over existing 
methods. Realistic finite element models can be introduced, without first solving large eigen- 
value problems to obtain the normal modes. This modeling can be done at the element level, 
using existing finite element libraries. No assembly of mass and stiffness matrices into global 
(system) matrices is required, hence there is no need to introduce special procedures for deal- 
ing with sparse matrices. Results indicate that the new computational scheme is capable of 
reproducing the energy exchange between the fluid and the structure with significantly less 
error than existing methods. 

Because no global mass and stiffness matrices are assembled for the structure, the actual 
dynamics simulated by the present method differs from what would be calculated by classical 
modal methods. The present method simulates the local force and momentum transfer 
between individual elements, in the time domain, and is therefore able to model wave propaga- 
tion phenomena that most modal methods are ill-equipped to do. Because the flutter mode 
emerges as part of the solution, as it would in a flutter test, there is no need to worry about 
"how many" and "which” modes should be kept in the flutter analysis. 

The importance of camber bending in certain transonic flutter problems has been illustrated, 
including a new trailing-edge "buzz” phenomenon that is predicted to occur if the chordwise 
stiffness is sufficiently low, and a nonclassical (inviscid) form of aileron buzz. Nonlinear 
aeroelastic modes, where two or more degrees-of-freedom exchange energy in a nonstationary 
manner during flutter, have been observed in numerical simulations of transonic flutter in cas- 
cades. No linear theory can be expected to provide an adequate theoretical basis for analyzing 
such phenomena. In some of these problems, many degrees of freedom participate, making 
classical modal method less efficient It is precisely in such problems that the advantages of the 
present method are brought out. For example, the possibility of considering traveling wave 
flutter modes without special advance set-up is an obvious advantage, as demonstrated in the 
transonic panel flutter solutions. 
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Fig. 2a Hybrid cascade mesh (reference channel, near field) 
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Fig. 2b Unstructured mesh used in panel flutter calculations (near field) 
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c) 


Fig. 3 a) Modified typical section wing model with camber bending 

b) Basic finite element 

c) Aileron hinge element 
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Fig. 5 Bifurcation diagram for the NACA 64A006 wing model, as calculated by present 
method and compared with Ref. 22: 

M « 0.87: A --- Ref. 22; Present method: A — t*/l = 0.05; O — t,A = °- 25 
M-0.92: ■ --- Ref. 22; Present method: □ — t,A=0.05; B ~t,A = 0.10 
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Fig. 6 Effect of chordwise flexibility on transonic limit cycle flutter at Mach 0.87 
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Fig. 10 Motion of supersonic pockets and shocks during symmetry-breaking bifurca- 
tion that gives rise to weak divergence and flutter-divergence interactions at Mach 0.80 
and p=75 
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Fig. 1 1 Strong flutter of NACA 0012 model at Mach 0.94 and U = 6.5, as calculated by 
present scheme. Classical scheme fails to predict flutter for this Mach number (see 
Fig. 12) and predicts divergence for U 2 8.3. Present scheme predicts no divergence- 
flutter interactions until U £ 1 1 



Fig. 12 Stable response of NACA 0012 model at Mach 0.94 and U = 8.0, as calculated 
by classical scheme. Classical scheme fails to predict flutter for this Mach number, up 
to the divergence boundary, in sharp contrast to the predictions of the new scheme 
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Fig. 14 Blade pressure coefficient C p during limit cycle flutter shown in Fig. 13 
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Fig 15 Nonlinear (nonclassical) bending-torsion flutter ot NACA 0006 cascade with 
6 = 45 deg; s 2 /c = 1 ; k, = 3.0; k„ = 0.375; EA at 0.512c 




Fig. 16 Highly nonlinear flutter at Mach 0.90 of cascade in Fig. 15. Note stabilizing 
effect of increasing Mach number from 0.85 to 0.90 
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Fig 17 Explosive flutter of a rotating NACA 0003 blade in an unstaggered cascade: a) 
plunging motion (at EA); b) angle of twist; c) Hamiltonian and work done by aero- 
dynamic forces (s/c = 1, = 0.85, o= 180 deg, L/C = 1) 
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Fig. 18 Transonic limit cycle panel flutter at Mach 1.0, illustrating highly irregular (but 
periodic) aeroelastic response and evidence of a traveling wave. Simply supported 
aluminum panel at 20,000 ft altitude; h/c = 0.004 
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Fig. 19 Transonic panel flutter at Mach 1.2 of same panel as in Fig. 18. Initial transient 
response is not shown 
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Fig. 20 Snapshot of panel deflection and corresponding Mach number distribution dur- 
ing limit cycle flutter at Mach 1 (Fig. 18). Flutter mode is a traveling wave and a strong 
shock is present on the panel surface 
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Fig. 21 Snapshots of panel flutter mode shape at different times during one cycle of 
flutter of panel at Mach 1 .2, corresponding to Fig. 19 


O. Bendiksen 35 



